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FOREWORD 


Stochastic  processes  are  utilized  to  represent  many  diverse 
phenomena.  This  report  describes  a  stochastic  process  generator 
useful  for  simulation  studies  with  the  help  of  the  computer. 

This  report  has  been  reviewed  by  John  H.  Walker,  Jr. 
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ABSTRACT 


A  method  has  been  derived  to  simulate  a  one -dimensional  stationary 
stochastic  process  with  a  given  autocorrelation  function  by  a  finite 
trigonometric  sum.  The  coefficients  of  the  latter  are  uncorrelated 
random  numbers.  A  rigorous  estimate  of  the  degree  of  approximation 
to  the  autocorrelation  function  is  given.  The  method  is  quite  general 
and  does  not  require  the  power  spectrum  to  be  rational. 
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1.  INTRODUCTION 


Stochastic  processes  (see  Refs.  1,  2)  are  used  to  represent 
phenomena  in  many  diverse  fields. 

It  is  of  importance  to  possess  numerical  algorithms  that 
generate  stochastic  processes  in  order  to  conduct  simulation 
studies  on  a  computer. 

In  this  report  we  adhere  to  the  convention  of  underlining  all 
random  quantities . 

The  method  to  be  outlined  uses  a  sequence  of  independent 
random  variables  to  generate  a  stationary  random  process  with  a 
specified  autocorrelation.  Thus  we  may  simulate  the 

process  on  a  computer  with  the  aid  of  a  "random  number  generator". 
This  report  confines  itself  to  discussing  one -dimensional  stochastic 
processes . 

2 .  Simulation  of  a  one -dimensional  weakly  stationary  stochastic 

process  with  given  autocorrelation. 

Let  the  mean  of  a  stationary  process  x(t)  be  zero  and  the 
autocorrelation  be  Rx(r).  The  power  spectrum  (see  ref.  1, 
p.  338)  is  given  by 
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The  inverse  relationship  is 

oo 

MT>  -  JT  /  el“T  s*(w)d"  (2) 
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We  have  here  and  forthwith  assumed  that  there  is  no  line  spectrum. 

We  may  now  construct  a  stochastic  process  x(t)  which  has  the 
given  mean,  possesses  approximately  the  given  autocorrelation  func 
tion  and  depends  only  on  a  discrete  number  of  random  variables. 

We  proceed  by  dividing  the  frequency  range  into  discrete  intervals 
Let  us  denote  the  size  of  the  interval  by  o>0.  This  increment  is 
to  be  small  in  some  sense  to  be  specified  later.  Define  the 
quantities 

nu0  +  i  "o 

«„  -  JW  (  Sx(6»dU  O) 

nw  -  ^  U> 

o  2  o 
n  =  0 ,  ±1 ,  ... 

These  quantities  are  obviously  non-negative  and  represent  the 
power  in  the  frequency  range  of  nw0  -  j  0)Q  to  na>Q  +  j  o>Q. 

Now  take  a  discrete  set  of  complex  random  variables  which  satisfy 
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Now,  of  course,  these  conditions  do  not  determine  the  random 
variables  £n  uniquely.  However  all  the  stochastic  processes  defined 


«<*)  -  2  £n  ein<Jot 

n  =  -oo 


possess  certain  common  properties.  We  have 
E  (x(t)  }  =  0  =  E  {  x(t) } 


E{x(t)2}  =  £  E  {  I  c  n  |  2  >  =  £  otn 
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Thus  the  processes  x(t)  have  the  same  mean  and  total  power  as  the 
process  x(t)  . 

Now  let  us  investigate  the  autocorrelation  of  x(t) .  We  have 
from  Eqs .7,8 

R£(t>  -  E  {x(t  +  T)x(t)) 
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Thus  the  x(t)  are  stationary  in  the  weak  sense,  but  their  auto¬ 
correlation  differs  from  R^T).  Let  us  .estimate  the  difference. 
How  well  does  R£(T)  approximate  ^(t)?  We  may  write 
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Each  term  in  Eq.  13  may  be  written 
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Therefore  we  may  write  for  the  difference  of  the  desired  auto¬ 
correlation  and  the  simulated  one 
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The  magnitude  of  the  final  expression  may  be  estimated  without 
difficulty.  When 

nWo  -  \  nuo  +  \  0)o  (16) 

we  have  the  inequality 
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Consequently  by  the  first  mean  value  theorem  of  the  integral 
calculus,  we  have  the  following  inequality 
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Equation  15  thus  yields  the  inequality, 

f°° 

rx(t)-ra(T)  ^  |  u>0\r\  ±-  J^sx(w)dw 

=  “  WQ|  Tl  1^(0) 

=  I  T  |  X  Total  Power  of  x(t)  (19) 


We  see  from  this  inequality  that  in  order  to  approximate  R  (t) 

x 

over  a  range  |T  |  ^  T,  we  must  choose  a  frequency  increment 

such  that 


(20) 


Furthermore  we  see  that  however  small  we  choose  a>o,  Ra(t)  will 
not  approximate  RX(T)  for  all  T.  This  is  connected  with  the 
fact  the  process  x(t)  involves  only  a  discrete  set  of  frequencies. 
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The  inequality  19  gives  us  the  effect  of  the  frequency 
discretization  on  the  autocorrelation  function.  For  the  purpose 
of  numerical  simulation  we  must  also  estimate  the  effect  of  trun¬ 
cating  the  series  in  Eq.  9.  In  contrast  to  the  frequency 
discretization,  the  truncation  affects  the  total  power  according 
to  Eq.  11.  If  we  truncate  after  |n|  =  N-l,.the  power  will  be 

reduced  according  to  Eq.  11  by 
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3 .  An  example 

Suppose  we  wish  to  simulate  a  stochastic  process  with  auto¬ 
correlation  function 

Rx(T)  =  <r2e"P|r|  (22) 


Its  power  spectrum  is  (see  Ref.  1,  p.  340) 


2  P(?2 
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The  total  power  is 
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We  may  define  a  correlation  distance  by 
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(25) 


Suppose  we  wish  to  represent  the  autocorrelation  function 
to  within  1006  %  of  the  total  power  over  a  distance  of  m 
correlation  distances.  The  maximum  lag  is  given  by 


m 


max  p 


(26) 


The  error  in  the  autocorrelation  function  due  to  discretization 
satisfies  according  to  Eq.  19 


^  m 
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The  error  due  to  truncation  satisfies  according  to  Eq.  21 


(28) 


(29) 
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and  therefore 


4.  Does  the  process  x(t)  approximate  the  process  x(t)? 

We  might  be  tempted  to  believe  that  a  process  &(t)  constructed 
in  the  aforementioned  way,  approximates  in  some  sense  the  process 
x(t) .  However  that  is  only  true  if 

r°°  sin(  i  (j  t) 

cn  =  j  x(t)  -  e"in u)Qt  dt  (37) 

-80  7T  t 

(see  Ref.  1,  p.  461).  These _cn  satisfy  all  the  conditions  expressed 
by  Eqs.  4-8  and  the  process  x(t)  then  approximates  x(t)  in  a  mean 
square  sense  as  discussed  in  the  reference. 

5.  Strict  sense  stationarity  and  normality 

It  may  be  useful  for  some  applications  to  be  able  to  generate 
stochastic  processes  that  are  stationary  in  the  strong  sense.  The 
stochastic  processes  x(t)  as  defined  by  Eqs.  3-9  are  in  general 
only  stationary  in  the  weak  sense.  We  may  show  this  by  regarding 
the  following  example:  Write 


c  = 
n 

i£n 

2n  e 

(38) 

■*n  = 

(39) 

in- 

-i-  n 

(40) 

Let  the  _pn  and  6n  be  mutually  independent  random  variables. 
Let  the  £n  have  the  probability  density  function 


f„  (.6)  =  —  (1  +  cos  30  ) 
^  2  IT 
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Let  the  _£n  satisfy 
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=  0 

(42) 
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n  =  1 ,  ... 
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n  =  0 ,  1,  ... 

(44) 
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(45) 

where  the  j3n  are  any  sequence  of  numbers  not  all  zero.  It 
follows  that  the  Eqs .  4-8  are  satisfied  and  the  process 

4(0-  £  £„einUot  (46) 

n  =  -oo 

is  weakly  stationary.  It  is  however  not  stationary  in  the  strong 
sense  as  may  be  shown  by  observing  that 


*(t)3 


E 


3in  U>0t 


(47) 


which  is  not  constant. 

It  may  be  shown  that  if  the  and  _£n  are  mutually 
independent,  the  satisfying  Eqs.  39,  42-44  and  the  ^  are 
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uniformly  distributed  in  -7 T  to  7T  and  satisfying  Eq.  40  the 
process  x(t)  of  Eq.  46  is  stationary  in  the  strict  sense. 

It  is  possible  to  generate  a  normal  process  by  the  method 
described  (for  the  existance  of  normal  processes  see  Ref.  2,  p.  72). 
In  this  case  weak  stationarity  implies  strong  s tationarity . 
Normality  is  often  assumed  if  the  process  to  be  simulated  is 
imperfectly  known.  Various  theoretical  results  concerning  sample 
distributions  may  be  derived  for  normal  process. 

To  generate  a  normal  process  let  us  write 


£n  =  |  (an  "  ibji) 

(48) 

in  £-n  3  ]ln 

=  -  b_n 

(49) 

We  may  rewrite  Eq.  9  as 

A  1  °° 

x(t)  =  j  2.0  +  E  (an  cos  nwQt  +  bn  sin  nUct) 
n  =  1 

(50) 

Let  us  choose  the  a_n  and  bn  as 

mutually  independent  normal 

random 

variables  satisfying 

E'{ib  }  =  0 

(51) 

E  {  an  }  =  E  {  bn 

}  =  0,  n  =  1,  ... 

(52) 

E(in}=  4«n> 

n  =  0 ,  1 ,  ... 

(53) 

E{bn2  }=  4oen, 

n  =  1 ,  ... 

(54) 
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These  equations  imply  Eqs .  4-8.  Now  since  x(t)  is  a  linear 
combination  of  normal  random  variables,  it  is  a  normal  process 
and  stationary  in  the  strict  sense. 

We  may  prove  that  the  normal  process  we  have  just  described 
has  the  property  of  "random  phase",  namely  the  random  variables 
defined  by 

tg£n  =  ^  (55) 


-7T  ^  $n  <  7 T  (56) 

are  mutually  independent  and  uniformly  distributed.  Indeed  the 
random  variable 

z  =  —  (57) 

£n 

is  distributed  according  to  Student* s  distribution  of  one  degree 
of  freedom  (Ref.  3,  p.  237). 

Student's  distribution  reduces  for  one  degree  of  freedom  to 
Cauchy's  distribution  (Ref.  3,  p.  246).  The  probability  density 
function  is 


f  (z)  =  —  - ± — - 

7T  i  +  z2 


(58) 


From  this  we  easily  derive  that  the  probability  density  function 


- 


(59) 


7f  ,  -7T  <  0  <  7T 

^  0 ,  otherwise 

Random  phases  are  often  made  a  requirement  in  representing  various 
phenomena  in  the  time  domain. 

6.  CONCLUSION 

We  have  presented  a  method  to  simulate  stationary  stochastic 
processes  numerically.  The  method  is  quite  general  and  does  not 
depend  on  the  power  spectrum  being  rational.  It  only  requires  a 
random  number  generator  and  adequate  computing  power. 
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